C     Last change:  ERB  10 Jul 2002    3:31 pm
      SUBROUTINE GWF1LPF1AL(ISUM,LCHK,LCVKA,LCSC1,LCSC2,LCHANI,
     1  LCVKCB,IN,NCOL,NROW,NLAY,IOUT,ILPFCB,LCWETD,
     2  HDRY,NPLPF,NCNFBD,LCLAYF,IREWND,ISUMI,LAYHDT,ITRSS,LCSV,ISEN)
C
C-----VERSION 11JAN2000 GWF1LPF1AL
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR LAYER PROPERTY FLOW PACKAGE
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------

      use altparam                               !jd

      INTEGER LAYHDT(NLAY)
      CHARACTER*14 LAYPRN(5),AVGNAM(3),TYPNAM(2),VKANAM(2),WETNAM(2),
     1            HANNAM
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA AVGNAM/'      HARMONIC','   LOGARITHMIC','     LOG-ARITH'/
      DATA TYPNAM/'      CONFINED','   CONVERTIBLE'/
      DATA VKANAM/'    VERTICAL K','    ANISOTROPY'/
      DATA WETNAM/'  NON-WETTABLE','      WETTABLE'/
      DATA HANNAM/'      VARIABLE'/
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
      IREWND=0
      ZERO=0.
C
C1------IDENTIFY PACKAGE
      WRITE(IOUT,1) IN
    1 FORMAT(1X,/1X,'LPF1 -- LAYER PROPERTY FLOW PACKAGE, VERSION 1',
     1', 1/11/2000',/,9X,'INPUT READ FROM UNIT ',I4)
C
C2------READ FIRST RECORD AND WRITE
      CALL URDCOM(IN,IOUT,LINE)
      LLOC=1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ILPFCB,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,3,I,HDRY,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NPLPF,R,IOUT,IN)
      IF(ILPFCB.LT.0) WRITE(IOUT,8)
    8 FORMAT(1X,'CONSTANT-HEAD CELL-BY-CELL FLOWS WILL BE PRINTED',
     1  ' WHEN ICLPFL IS NOT 0')
      IF(ILPFCB.GT.0) WRITE(IOUT,9) ILPFCB
    9 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
      if(hdrybot.eq.0)then                                             !jd
      WRITE(IOUT,11) HDRY
   11 FORMAT(1X,'HEAD AT CELLS THAT CONVERT TO DRY=',1PG13.5)
      else                                                             !jd
        write(iout,1111)                                               !jd
 1111   FORMAT(1X,'HEAD AT CELLS THAT CONVERT TO DRY= layer bottom ',  !jd
     +'elevation.')                                                    !jd
      end if                                                           !jd
C
C  Allocate space in parameter-information arrays
      CALL UPARARRAL(0,IOUT,LINE,NPLPF)
C
C3------STOP THE SIMULATION IF THERE ARE MORE THAN 200 LAYERS.
   39 IF(NLAY.GT.200) THEN
         WRITE(IOUT,41)
   41    FORMAT(1X,/1X,'YOU HAVE SPECIFIED MORE THAN 200 MODEL LAYERS'/
     1 1X,'SPACE IS RESERVED FOR A MAXIMUM OF 200 LAYERS IN LPF ARRAYS')
        call stopfile ! emrl
         STOP
      END IF
C
C4------READ LAYTYP, LAYAVG, CHANI, LAYVKA, LAYWET.
      READ(IN,*) (LAYTYP(K),K=1,NLAY)
      READ(IN,*) (LAYAVG(K),K=1,NLAY)
      READ(IN,*) (CHANI(K),K=1,NLAY)
      READ(IN,*) (LAYVKA(K),K=1,NLAY)
      READ(IN,*) (LAYWET(K),K=1,NLAY)
C
C4A-----PRINT A TABLE OF VALUES FOR LAYTYP, LAYAVG, CHANI, LAYVKA, LAYWET.
      WRITE(IOUT,47)
   47 FORMAT(1X,/3X,'LAYER FLAGS:',/1X,
     1 'LAYER       LAYTYP        LAYAVG         CHANI ',
     2 '       LAYVKA        LAYWET',/1X,75('-'))
      DO 50 K=1,NLAY
      WRITE(IOUT,48) K,LAYTYP(K),LAYAVG(K),CHANI(K),LAYVKA(K),LAYWET(K)
   48 FORMAT(1X,I4,2I14,1PE14.3,2I14)
C     SET GLOBAL HEAD-DEPENDENT THICKNESS INDICATOR
      IF (LAYTYP(K).NE.0) THEN
        LAYHDT(K)=1
      ELSE
        LAYHDT(K)=0
      ENDIF
   50 CONTINUE
C
C4B-----BASED ON LAYTYP, LAYAVG, CHANI, LAYWET, COUNT THE NUMBER OF EACH
C4B-----TYPE OF 2-D ARRAY; CHECK VALUES FOR CONSISTENCY; AND SETUP
C4B-----POINTERS IN LAYTYP, CHANI, AND LAYWET FOR CONVENIENT ACCESS
C4B-----TO SC2, HANI, and WETDRY.  PRINT INTERPRETED VALUES OF FLAGS.
      NCNVRT=0
      NHANI=0
      NWETD=0
      WRITE(IOUT,57)
   57 FORMAT(1X,/3X,'INTERPRETATION OF LAYER FLAGS:',/1X,
     1  '                       INTERBLOCK     HORIZONTAL',
     2  '    DATA IN',/1X,
     3  '        LAYER TYPE   TRANSMISSIVITY   ANISOTROPY',
     4  '   ARRAY VKA   WETTABILITY',/1X,
     5  'LAYER      (LAYTYP)      (LAYAVG)       (CHANI)',
     6  '      (LAYVKA)      (LAYWET)',/1X,75('-'))
      DO 100 K=1,NLAY
      IF(LAYTYP(K).NE.0) THEN
         NCNVRT=NCNVRT+1
         LAYTYP(K)=NCNVRT
      END IF
      IF(CHANI(K).LE.ZERO) THEN
         NHANI=NHANI+1
         CHANI(K)=-NHANI
      END IF
      IF(LAYWET(K).NE.0) THEN
         IF(LAYTYP(K).EQ.0) THEN
            WRITE(IOUT,*)
     1          ' LAYWET is not 0 and LAYTYP is 0 for layer:',K
            WRITE(IOUT,*) ' LAYWET must be 0 if LAYTYP is 0'
        call stopfile ! emrl
            STOP
         ELSE
            NWETD=NWETD+1
            LAYWET(K)=NWETD
         END IF
      END IF
      IF(LAYAVG(K).LT.0 .OR. LAYAVG(K).GT.2) THEN
         WRITE(IOUT,74) LAYAVG(K)
   74    FORMAT(1X,I8,
     1    ' IS AN INVALID LAYAVG VALUE -- MUST BE 0, 1, or 2')
        call stopfile ! emrl
         STOP
      END IF
      LAYPRN(1)=TYPNAM(1)
      IF(LAYTYP(K).NE.0) LAYPRN(1)=TYPNAM(2)
      LAYPRN(2)=AVGNAM(LAYAVG(K)+1)
      IF(CHANI(K).LE.0) THEN
         LAYPRN(3)=HANNAM
      ELSE
         WRITE(LAYPRN(3),'(1PE14.3)') CHANI(K)
      END IF
      LAYPRN(4)=VKANAM(1)
      IF(LAYVKA(K).NE.0) LAYPRN(4)=VKANAM(2)
      LAYPRN(5)=WETNAM(1)
      IF(LAYWET(K).NE.0) LAYPRN(5)=WETNAM(2)
      WRITE(IOUT,78) K,(LAYPRN(I),I=1,5)
   78 FORMAT(1X,I4,5A)
  100 CONTINUE
C
C5------COMPUTE THE NUMBER OF CELLS IN THE ENTIRE GRID AND IN ONE LAYER.
      NRC=NROW*NCOL
      ISIZ=NRC*NLAY
C
C6------ALLOCATE SPACE FOR ARRAYS.
      ISOLD=ISUM
      ISOLDI=ISUMI
      LCVKA=ISUM
      ISUM=ISUM+ISIZ
      LCVKCB=ISUM
      ISUM=ISUM+NRC*NCNFBD
      LCSV = ISUM
      IF(ISEN.GT.0) ISUM=ISUM+NRC*(NLAY-1)
      LCSC1=ISUM
      IF(ITRSS.NE.0) ISUM=ISUM+ISIZ
      LCSC2=ISUM
      IF(ITRSS.NE.0) ISUM=ISUM+NRC*NCNVRT
      LCHANI=ISUM
      ISUM=ISUM+NRC*NHANI
      LCWETD=ISUM
      ISUM=ISUM+NRC*NWETD
      LCLAYF=ISUMI
      ISUMI=ISUMI+NLAY*6
      LCHK=ISUM
      ISUM=ISUM+ISIZ
C
C7------PRINT THE AMOUNT OF SPACE USED BY THE LPF PACKAGE.
      ISP=ISUM-ISOLD
      WRITE(IOUT,101) ISP
  101 FORMAT(1X,/1X,I10,' ELEMENTS IN X ARRAY ARE USED BY LPF')
      ISP=ISUMI-ISOLDI
      WRITE(IOUT,102) ISP
  102 FORMAT(1X,I10,' ELEMENTS IN IX ARRAY ARE USED BY LPF')
C
C8------RETURN.
      RETURN
      END
      SUBROUTINE GWF1LPF1RQ(HK,VKA,
     1 VKCB,HANI,SC1,SC2,IN,ITRSS,NCOL,NROW,NLAY,IOUT,WETDRY,NPLPF,
     2 WETFCT,IWETIT,IHDWET,LAYFLG,BOTM,NBOTM,DELR,DELC,ITERP,INAMLOC)
C
C-----VERSION 11JAN2000 GWF1LPF1RQ
C     ******************************************************************
C     READ DATA FOR LAYER PROPERTY FLOW PACKAGE
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      INCLUDE 'param.inc'
      CHARACTER*24 ANAME(8)
      CHARACTER*4 PTYP
C
      DIMENSION HK(NCOL,NROW,NLAY),LAYFLG(6,NLAY),
     1    VKA(NCOL,NROW,NLAY),VKCB(NCOL,NROW,NLAY),
     2    HANI(NCOL,NROW,NLAY),SC1(NCOL,NROW,NLAY),
     3    SC2(NCOL,NROW,NLAY),WETDRY(NCOL,NROW,NLAY),
     4    BOTM(NCOL,NROW,0:NBOTM),DELR(NCOL),DELC(NROW)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA ANAME(1) /'   HYD. COND. ALONG ROWS'/
      DATA ANAME(2) /'  HORIZ. ANI. (COL./ROW)'/
      DATA ANAME(3) /'     VERTICAL HYD. COND.'/
      DATA ANAME(4) /' HORIZ. TO VERTICAL ANI.'/
      DATA ANAME(5) /'QUASI3D VERT. HYD. COND.'/
      DATA ANAME(6) /'        SPECIFIC STORAGE'/
      DATA ANAME(7) /'          SPECIFIC YIELD'/
      DATA ANAME(8) /'        WETDRY PARAMETER'/
C     ------------------------------------------------------------------
  500 FORMAT(/,
     &' ERROR: WHEN A HANI PARAMETER IS USED, CHANI FOR ALL LAYERS',/,
     &' MUST BE LESS THAN OR EQUAL TO 0.0 -- STOP EXECUTION',
     &' (GWF1LPF1RQ)')
C
      ZERO=0.
C
      IWDFLG=0
      DO 10 K=1,NLAY
      IF(LAYWET(K).NE.0) IWDFLG=IWDFLG+1
   10 CONTINUE
      IF(IWDFLG.EQ.0) THEN
         WRITE(IOUT,11)
   11    FORMAT(1X,/,1X,'WETTING CAPABILITY IS NOT ACTIVE IN ANY LAYER')
      ELSE
         WRITE(IOUT,12) IWDFLG
   12    FORMAT(1X,/,1X,'WETTING CAPABILITY IS ACTIVE IN',I4,' LAYERS')
         IWDFLG=1
         READ(IN,*) WETFCT,IWETIT,IHDWET
         IF(IWETIT.LE.0) IWETIT=1
         WRITE(IOUT,*) ' WETTING FACTOR=',WETFCT
         WRITE(IOUT,*) ' WETTING ITERATION INTERVAL=',IWETIT
         WRITE(IOUT,*) ' IHDWET=',IHDWET
      END IF
C
C-------READ NAMED PARAMETERS
      NPHK=0
      NPVKCB=0
      NPVK=0
      NPVANI=0
      NPSS=0
      NPSY=0
c%ERA%HANI
      NPHANI=0
      IF(NPLPF.GT.0) THEN
         WRITE(IOUT,15)
   15    FORMAT(/,' PARAMETERS DEFINED IN THE LPF PACKAGE')
         DO 20 K=1,NPLPF
         CALL UPARARRRP(IN,IOUT,N,1,PTYP,ITERP,0,INAMLOC)
C   Note that NPHK and the other NP variables in
C   this group are used only as flags, not counts
         IF(PTYP.EQ.'HK') THEN
            NPHK=1
         ELSE IF(PTYP.EQ.'HANI') THEN
C           ENSURE THAT ALL CHANI <= 0
            DO 17 I = 1, NLAY
              IF (CHANI(I).GT.0.0) THEN
                WRITE(IOUT,500)
        call stopfile ! emrl
                STOP
              ENDIF
   17       CONTINUE
            NPHANI=1
         ELSE IF(PTYP.EQ.'VKCB') THEN 
            NPVKCB=1
         ELSE IF(PTYP.EQ.'VK') THEN
            NPVK=1
            CALL SGWF1LPF1CK(IOUT,N,'VK  ')
         ELSE IF(PTYP.EQ.'VANI') THEN
            NPVANI=1
            CALL SGWF1LPF1CK(IOUT,N,'VANI')
         ELSE IF(PTYP.EQ.'SS') THEN
            NPSS=1
         ELSE IF(PTYP.EQ.'SY') THEN
            NPSY=1
         ELSE
            WRITE(IOUT,*) ' Invalid parameter type for LPF Package'
        call stopfile ! emrl
            STOP
         END IF
C  Make the parameter active for all stress periods
         IACTIVE(N)=-1
   20    CONTINUE
      END IF
C
C2------DEFINE DATA FOR EACH LAYER -- VIA READING OR NAMED PARAMETERS.
      DO 200 K=1,NLAY
      KK=K
C
C2B-----DEFINE HORIZONTAL HYDRAULIC CONDUCTIVITY (HK)
      IF(NPHK.EQ.0) THEN
         CALL U2DREL(HK(1,1,KK),ANAME(1),NROW,NCOL,KK,IN,IOUT)
      ELSE
         READ(IN,*) LAYFLG(1,K)
         WRITE(IOUT,101) ANAME(1),K,LAYFLG(1,K)
  101    FORMAT(1X,/1X,A,' FOR LAYER',I4,
     1   ' WILL BE DEFINED BY PARAMETERS',/1X,'(PRINT FLAG=',I4,')')
      END IF
C
C2C-----READ HORIZONTAL ANISOTROPY IF CHANI IS NON-ZERO
      IF(CHANI(K).LE.ZERO) THEN
         KHANI=-CHANI(K)
C%ERA%HANI
        IF(NPHANI.EQ.0) THEN
           CALL U2DREL(HANI(1,1,KHANI),ANAME(2),NROW,NCOL,KK,IN,IOUT)
        ELSE
           READ(IN,*) LAYFLG(6,K)
           WRITE(IOUT,101) ANAME(2),K,LAYFLG(6,K)
        END IF
      END IF
C
C2D-----DEFINE VERTICAL HYDRAULIC CONDUCTIVITY OR HORIZONTAL TO VERTICAL
C2D-----ANISOTROPY (VKA).
      IANAME=3
      IF(LAYVKA(K).NE.0) IANAME=4
      IF(NPVK.EQ.0 .AND. NPVANI.EQ.0) THEN
         CALL U2DREL(VKA(1,1,KK),ANAME(IANAME),NROW,NCOL,KK,IN,IOUT)
      ELSE
         READ(IN,*) LAYFLG(2,K)
         WRITE(IOUT,101) ANAME(IANAME),K,LAYFLG(2,K)
      END IF
C
C2E-----DEFINE SPECIFIC STORAGE IN ARRAY SC1 IF TRANSIENT.
      IF(ITRSS.NE.0) THEN
         IF(NPSS.EQ.0) THEN
            CALL U2DREL(SC1(1,1,KK),ANAME(6),NROW,NCOL,KK,IN,IOUT)
            CALL SGWF1LPF1SC(SC1(1,1,KK),NCOL,NROW,KK,BOTM,NBOTM,
     1            DELR,DELC,1)
         ELSE
            READ(IN,*) LAYFLG(3,K)
            WRITE(IOUT,101) ANAME(6),K,LAYFLG(3,K)
         END IF
      END IF
C
C2F-----DEFINE SPECIFIC YIELD IN ARRAY SC2 IF TRANSIENT AND LAYER IS
C2F-----IS CONVERTIBLE.
      IF(LAYTYP(K).NE.0) THEN
         IF(ITRSS.NE.0) THEN
            IF(NPSY.EQ.0) THEN
               CALL U2DREL(SC2(1,1,LAYTYP(K)),ANAME(7),NROW,NCOL,KK,IN,
     1                 IOUT)
               CALL SGWF1LPF1SC(SC2(1,1,LAYTYP(K)),NCOL,NROW,KK,BOTM,
     1                 NBOTM,DELR,DELC,0)
            ELSE
               READ(IN,*) LAYFLG(4,K)
               WRITE(IOUT,101) ANAME(7),K,LAYFLG(4,K)
            END IF
         END IF
      END IF
C
C2G-----READ CONFINING BED VERTICAL HYDRAULIC CONDUCTIVITY (VKCB)
      IF(LAYCBD(K).NE.0) THEN
         IF(NPVKCB.EQ.0) THEN
            CALL U2DREL(VKCB(1,1,LAYCBD(K)),ANAME(5),NROW,NCOL,KK,IN,
     1             IOUT)
         ELSE
            READ(IN,*) LAYFLG(5,K)
            WRITE(IOUT,101) ANAME(5),K,LAYFLG(5,K)
         END IF
      END IF
C
C2H-----READ WETDRY CODES IF WETTING CAPABILITY HAS BEEN INVOKED
C2H-----(LAYWET NOT 0).
      IF(LAYWET(K).NE.0) THEN
         CALL U2DREL(WETDRY(1,1,LAYWET(K)),ANAME(8),NROW,NCOL,KK,IN,
     1            IOUT)
      END IF
  200 CONTINUE
C
C4------RETURN
      RETURN
      END
      SUBROUTINE GWF1LPF1SP(IBOUND,HNEW,CR,CC,CV,DELR,DELC,BOTM,HK,VKA,
     1 VKCB,HANI,SC1,SC2,ITRSS,NCOL,NROW,NLAY,IOUT,WETDRY,NPLPF,
     2 NBOTM,RMLT,IZON,NMLTAR,NZONAR,LAYFLG,BUFF,ITERPK)
C
C-----VERSION 11JAN2000 GWF1LPF1SP
C     ******************************************************************
C     SUBSTITUTE AND PREPARE DATA FOR LAYER PROPERTY FLOW PACKAGE
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*24 ANAME(8)
      DOUBLE PRECISION HNEW
C
      DIMENSION IBOUND(NCOL,NROW,NLAY),HNEW(NCOL,NROW,NLAY),
     1    CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     2    CV(NCOL,NROW,NLAY),DELR(NCOL),DELC(NROW),
     3    BOTM(NCOL,NROW,0:NBOTM),HK(NCOL,NROW,NLAY),
     4    VKA(NCOL,NROW,NLAY),VKCB(NCOL,NROW,NLAY),
     5    HANI(NCOL,NROW,NLAY),SC1(NCOL,NROW,NLAY),
     6    SC2(NCOL,NROW,NLAY),WETDRY(NCOL,NROW,NLAY),
     7    RMLT(NCOL,NROW,NMLTAR),IZON(NCOL,NROW,NZONAR),
     8    LAYFLG(6,NLAY),BUFF(NCOL,NROW)
C
      INCLUDE 'param.inc'
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA ANAME(1) /'   HYD. COND. ALONG ROWS'/
      DATA ANAME(2) /'  HORIZ. ANI. (COL./ROW)'/
      DATA ANAME(3) /'     VERTICAL HYD. COND.'/
      DATA ANAME(4) /' HORIZ. TO VERTICAL ANI.'/
      DATA ANAME(5) /'QUASI3D VERT. HYD. COND.'/
      DATA ANAME(6) /'        SPECIFIC STORAGE'/
      DATA ANAME(7) /'          SPECIFIC YIELD'/
      DATA ANAME(8) /'        WETDRY PARAMETER'/
C     ------------------------------------------------------------------
C
      ZERO=0.
C
C-------LOOK FOR NAMED PARAMETERS FOR LPF
      NPHK=0
      NPVKCB=0
      NPVK=0
      NPVANI=0
      NPSS=0
      NPSY=0
C%ERA%HANI
      NPHANI=0
      IF(NPLPF.GT.0) THEN
         DO 20 K=1,MXPAR
C   Note that NPHK and the other NP variables in
C   this group are used only as flags, not counts
         IF(PARTYP(K).EQ.'HK') THEN
            NPHK=1
         ELSE IF(PARTYP(K).EQ.'HANI') THEN
            NPHANI=1
         ELSE IF(PARTYP(K).EQ.'VKCB') THEN
            NPVKCB=1
         ELSE IF(PARTYP(K).EQ.'VK') THEN
            NPVK=1
         ELSE IF(PARTYP(K).EQ.'VANI') THEN
            NPVANI=1
         ELSE IF(PARTYP(K).EQ.'SS') THEN
            NPSS=1
         ELSE IF(PARTYP(K).EQ.'SY') THEN
            NPSY=1
         END IF
   20    CONTINUE
      END IF
C
C2------DEFINE DATA FOR NAMED PARAMETERS.
      DO 200 K=1,NLAY
      KK=K
C
C2B-----DEFINE HORIZONTAL HYDRAULIC CONDUCTIVITY (HK)
      IF(NPHK.NE.0) THEN
         CALL UPARARRSUB1(HK(1,1,KK),NCOL,NROW,KK,'HK',
     1      IOUT,ANAME(1),LAYFLG(1,KK),RMLT,IZON,NMLTAR,NZONAR)
         IF (ITERPK.EQ.1) CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,NCOL,
     &                                   NLAY,NROW,NZONAR,'HK  ')
      ENDIF
C%ERA%HANI
C
C2C-----DEFINE HORIZONTAL ANISOTROPY (HANI) IF CHANI IS NON-ZERO
      IF(CHANI(K).LE.ZERO.AND.NPHANI.NE.0) THEN
        KHANI=-CHANI(K)
        CALL UPARARRSUB1(HANI(1,1,KHANI),NCOL,NROW,KK,'HANI',
     &      IOUT,ANAME(1),LAYFLG(6,KK),RMLT,IZON,NMLTAR,NZONAR)
        IF (ITERPK.EQ.1) CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,NCOL,
     &                                  NLAY,NROW,NZONAR,'HANI')
      ENDIF
C
C2D-----DEFINE VERTICAL HYDRAULIC CONDUCTIVITY OR HORIZONTAL TO VERTICAL
C2D-----ANISOTROPY (VKA ARRAY).
      IF (NPVK.NE.0 .OR. NPVANI.NE.0) THEN
         IF (LAYVKA(K).EQ.0) THEN
            IANAME=3
            CALL UPARARRSUB1(VKA(1,1,KK),NCOL,NROW,KK,'VK',IOUT,
     &                       ANAME(IANAME),LAYFLG(2,KK),RMLT,IZON,
     &                       NMLTAR,NZONAR)
            IF (ITERPK.EQ.1) CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,
     &                                      NCOL,NLAY,NROW,NZONAR,
     &                                      'VK  ')
         ELSE
            IANAME=4
            CALL UPARARRSUB1(VKA(1,1,KK),NCOL,NROW,KK,'VANI',IOUT,
     &                       ANAME(IANAME),LAYFLG(2,KK),RMLT,IZON,
     &                       NMLTAR,NZONAR)
            IF (ITERPK.EQ.1) CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,
     &                                      NCOL,NLAY,NROW,NZONAR,
     &                                      'VANI')
         ENDIF
      ENDIF
C
C2E-----DEFINE SPECIFIC STORAGE IN ARRAY SC1 IF TRANSIENT.
      IF(ITRSS.NE.0) THEN
         IF(NPSS.NE.0) THEN
            CALL UPARARRSUB1(SC1(1,1,KK),NCOL,NROW,KK,'SS',
     1      IOUT,ANAME(6),LAYFLG(3,KK),RMLT,IZON,NMLTAR,NZONAR)
            CALL SGWF1LPF1SC(SC1(1,1,KK),NCOL,NROW,KK,BOTM,NBOTM,
     1            DELR,DELC,1)
            IF (ITERPK.EQ.1)
     &          CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,NCOL,NLAY,NROW,
     &                         NZONAR,'SS  ')
         END IF
      END IF
C
C2F-----DEFINE SPECIFIC YIELD IN ARRAY SC2 IF TRANSIENT AND LAYER IS
C2F-----IS CONVERTIBLE.
      IF(LAYTYP(K).NE.0) THEN
         IF(ITRSS.NE.0) THEN
            IF(NPSY.NE.0) THEN
               CALL UPARARRSUB1(SC2(1,1,LAYTYP(K)),NCOL,
     1         NROW,KK,'SY',IOUT,ANAME(7),LAYFLG(4,KK),RMLT,IZON,
     2         NMLTAR,NZONAR)
               CALL SGWF1LPF1SC(SC2(1,1,LAYTYP(K)),NCOL,NROW,KK,BOTM,
     1             NBOTM,DELR,DELC,0)
               IF (ITERPK.EQ.1)
     &             CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,NCOL,NLAY,
     &                            NROW,NZONAR,'SY  ')
            END IF
         END IF
      END IF
C
C2G-----DEFINE CONFINING BED VERTICAL HYDRAULIC CONDUCTIVITY (VKCB)
      IF(LAYCBD(K).NE.0) THEN
         IF(NPVKCB.NE.0) THEN
            CALL UPARARRSUB1(VKCB(1,1,LAYCBD(K)),NCOL,NROW,KK,
     1         'VKCB',IOUT,ANAME(5),LAYFLG(5,KK),RMLT,IZON,
     2         NMLTAR,NZONAR)
            IF (ITERPK.EQ.1)
     &          CALL UPARARRCK(BUFF,IBOUND,IOUT,IZON,K,NCOL,NLAY,NROW,
     &                         NZONAR,'VKCB')
         ENDIF
      END IF
C
  200 CONTINUE
C
C3------PREPARE AND CHECK LPF DATA.
      CALL SGWF1LPF1N(HNEW,IBOUND,CR,CC,CV,HK,HANI,VKA,VKCB,DELR,
     1       DELC,NCOL,NROW,NLAY,IOUT,WETDRY,BOTM,NBOTM)
C
C4------RETURN
      RETURN
      END
      SUBROUTINE GWF1LPF1AD(IBOUND,HOLD,BOTM,WETDRY,ISS,NCOL,NROW,NLAY,
     1            NBOTM)
C
C-----VERSION 11JAN2000 GWF1LPF1AD
C     ******************************************************************
C     SET HOLD TO BOTM WHENEVER A WETTABLE CELL IS DRY
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DIMENSION IBOUND(NCOL,NROW,NLAY),HOLD(NCOL,NROW,NLAY),
     1          BOTM(NCOL,NROW,0:NBOTM),WETDRY(NCOL,NROW,NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
C
C1------RETURN IF STEADY STATE.
      IF(ISS.NE.0) RETURN
C
C2------LOOP THROUGH ALL LAYERS TO SET HOLD=BOT IF A WETTABLE CELL IS DRY
      ZERO=0.
      DO 100 K=1,NLAY
C
C2A-----SKIP LAYERS THAT CANNOT CONVERT BETWEEN WET AND DRY
      IF(LAYWET(K).EQ.0) GO TO 100
      DO 90 I=1,NROW
      DO 90 J=1,NCOL
C
C2B-----SKIP CELLS THAT ARE CURRENTLY WET OR ARE NOT WETTABLE
      IF(IBOUND(J,I,K).NE.0) GO TO 90
      IF(WETDRY(J,I,LAYWET(K)).EQ.ZERO) GO TO 90
C
C2C-----SET HOLD=BOT
      HOLD(J,I,K)=BOTM(J,I,LBOTM(K))
   90 CONTINUE
  100 CONTINUE
C
C3-----RETURN
      RETURN
      END
      SUBROUTINE GWF1LPF1FM(HCOF,RHS,HOLD,SC1,HNEW,IBOUND,CR,CC,CV,HK,
     1      HANI,VKA,BOTM,SC2,DELR,DELC,DELT,ISS,KITER,KSTP,KPER,NCOL,
     2      NROW,NLAY,IOUT,WETDRY,WETFCT,IWETIT,IHDWET,HDRY,NBOTM,VKCB)
C-----VERSION 11JAN2000 GWF1LPF1FM
C     ******************************************************************
C     ADD LEAKAGE CORRECTION AND STORAGE TO HCOF AND RHS, AND CALCULATE
C     CONDUCTANCE AS REQUIRED.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW
C
      DIMENSION HCOF(NCOL,NROW,NLAY),RHS(NCOL,NROW,NLAY),
     1    HOLD(NCOL,NROW,NLAY),SC1(NCOL,NROW,NLAY),HNEW(NCOL,NROW,NLAY),
     2    IBOUND(NCOL,NROW,NLAY),CR(NCOL,NROW,NLAY),
     3    CC(NCOL,NROW,NLAY),CV(NCOL,NROW,NLAY),HK(NCOL,NROW,NLAY),
     4    HANI(NCOL,NROW,NLAY),VKA(NCOL,NROW,NLAY),
     4    BOTM(NCOL,NROW,0:NBOTM),SC2(NCOL,NROW,NLAY),DELR(NCOL),
     5    DELC(NROW),WETDRY(NCOL,NROW,NLAY),VKCB(NCOL,NROW,NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
      ONE=1.
C
C1------FOR EACH LAYER: IF CONVERTIBLE, CALCULATE CONDUCTANCES.
      DO 100 K=1,NLAY
      KK=K
      IF(LAYTYP(K).NE.0)
     1    CALL SGWF1LPF1HCOND(HNEW,IBOUND,CR,CC,HK,HANI,DELR,DELC,BOTM,
     2      NBOTM,KK,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT,WETDRY,
     3      WETFCT,IWETIT,IHDWET,HDRY)
  100 CONTINUE
      DO 101 K=1,NLAY
      KK=K
      IF(K.NE.NLAY) THEN
         IF(LAYTYP(K).NE.0 .OR. LAYTYP(K+1).NE.0)
     1      CALL SGWF1LPF1VCOND(CV,HK,VKA,VKCB,IBOUND,BOTM,NBOTM,KK,
     2         NCOL,NROW,NLAY,HNEW,DELR,DELC,IOUT)
      END IF
  101 CONTINUE
C
C2------IF THE STRESS PERIOD IS TRANSIENT, ADD STORAGE TO HCOF AND RHS
      IF(ISS.EQ.0) THEN
         TLED=ONE/DELT
         DO 200 K=1,NLAY
C
C3------SEE IF THIS LAYER IS CONVERTIBLE OR NON-CONVERTIBLE.
         IF(LAYTYP(K).EQ.0) THEN
C4------NON-CONVERTIBLE LAYER, SO USE PRIMARY STORAGE
            DO 140 I=1,NROW
            DO 140 J=1,NCOL
            IF(IBOUND(J,I,K).LE.0) GO TO 140
            RHO=SC1(J,I,K)*TLED
            HCOF(J,I,K)=HCOF(J,I,K)-RHO
            RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K)
  140       CONTINUE
         ELSE
C
C5------A CONVERTIBLE LAYER, SO CHECK OLD AND NEW HEADS TO DETERMINE
C5------WHEN TO USE PRIMARY AND SECONDARY STORAGE
            DO 180 I=1,NROW
            DO 180 J=1,NCOL
C
C5A-----IF THE CELL IS EXTERNAL THEN SKIP IT.
            IF(IBOUND(J,I,K).LE.0) GO TO 180
            TP=BOTM(J,I,LBOTM(K)-1)
            RHO2=SC2(J,I,LAYTYP(K))*TLED
            RHO1=SC1(J,I,K)*TLED
C
C5B-----FIND STORAGE FACTOR AT START OF TIME STEP.
            SOLD=RHO2
            IF(HOLD(J,I,K).GT.TP) SOLD=RHO1
C
C5C-----FIND STORAGE FACTOR AT END OF TIME STEP.
            HTMP=HNEW(J,I,K)
            SNEW=RHO2
            IF(HTMP.GT.TP) SNEW=RHO1
C
C5D-----ADD STORAGE TERMS TO RHS AND HCOF.
            HCOF(J,I,K)=HCOF(J,I,K)-SNEW
            RHS(J,I,K)=RHS(J,I,K) - SOLD*(HOLD(J,I,K)-TP) - SNEW*TP
C
  180       CONTINUE
         END IF
C
  200    CONTINUE
      END IF
C
C6------FOR EACH LAYER DETERMINE IF CORRECTION TERMS ARE NEEDED FOR
C6------FLOW DOWN INTO PARTIALLY SATURATED LAYERS.
      DO 300 K=1,NLAY
C
C7------SEE IF CORRECTION IS NEEDED FOR LEAKAGE FROM ABOVE.
      IF(LAYTYP(K).NE.0 .AND. K.NE.1) THEN
C
C7A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.
         DO 220 I=1,NROW
         DO 220 J=1,NCOL
C
C7B-----IF THE CELL IS EXTERNAL(IBOUND<=0) THEN SKIP IT.
         IF(IBOUND(J,I,K).LE.0) GO TO 220
         HTMP=HNEW(J,I,K)
C
C7C-----IF HEAD IS ABOVE TOP THEN CORRECTION NOT NEEDED
         TOP=BOTM(J,I,LBOTM(K)-1)
         IF(HTMP.GE.TOP) GO TO 220
C
C7D-----WITH HEAD BELOW TOP ADD CORRECTION TERMS TO RHS.
         RHS(J,I,K)=RHS(J,I,K) + CV(J,I,K-1)*(TOP-HTMP)
  220    CONTINUE
      END IF
C
C8------SEE IF THIS LAYER MAY NEED CORRECTION FOR LEAKAGE TO BELOW.
      IF(K.EQ.NLAY) GO TO 300
      IF(LAYTYP(K+1).NE.0) THEN
C
C8A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.
         DO 280 I=1,NROW
         DO 280 J=1,NCOL
C
C8B-----IF CELL IS EXTERNAL (IBOUND<=0) THEN SKIP IT.
         IF(IBOUND(J,I,K).LE.0) GO TO 280
C
C8C-----IF HEAD IN THE LOWER CELL IS LESS THAN TOP ADD CORRECTION
C8C-----TERM TO RHS.
         HTMP=HNEW(J,I,K+1)
         TOP=BOTM(J,I,LBOTM(K+1)-1)
         IF(HTMP.LT.TOP) RHS(J,I,K)=RHS(J,I,K)- CV(J,I,K)*(TOP-HTMP)
  280    CONTINUE
      END IF
C
  300 CONTINUE
C
C9------RETURN
      RETURN
      END
      SUBROUTINE SGWF1LPF1N(HNEW,IBOUND,CR,CC,CV,HK,HANI,VKA,VKCB,
     1        DELR,DELC,NCOL,NROW,NLAY,IOUT,WETDRY,BOTM,NBOTM)
C
C-----VERSION 13JAN2000 SGWF1LPF1N
C     ******************************************************************
C     INITIALIZE AND CHECK LPF DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DOUBLE PRECISION HNEW,HCNV
C
      DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     1    CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     2    CV(NCOL,NROW,NLAY),HK(NCOL,NROW,NLAY),
     3    HANI(NCOL,NROW,NLAY),VKA(NCOL,NROW,NLAY),VKCB(NCOL,NROW,NLAY),
     4    DELR(NCOL),DELC(NROW),WETDRY(NCOL,NROW,NLAY),
     5    BOTM(NCOL,NROW,0:NBOTM)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
      ZERO=0.
      HCNV=888.88
C
C1------INSURE THAT EACH ACTIVE CELL HAS AT LEAST ONE NON-ZERO
C1------TRANSMISSIVE PARAMETER.
      DO 60 K=1,NLAY
      IF(LAYWET(K).NE.0) THEN
C
C2------WETTING IS ACTIVE.
         DO 40 I=1,NROW
         DO 40 J=1,NCOL
         IF(IBOUND(J,I,K).EQ.0 .AND. WETDRY(J,I,LAYWET(K)).EQ.ZERO)
     1                GO TO 40
C
C2A-----CHECK HORIZONTAL HYDRAULIC CONDUCTIVITY (HK).
         IF(HK(J,I,K).NE.ZERO) GO TO 40
C
C2B-----CHECK VERTICAL HYDRAULIC CONDUCTIVITY AND CONFINING BED
C2B-----VERTICAL HYDRAULIC CONDUCTIVITY.
         IF(NLAY.GT.1) THEN
            IF(VKA(J,I,K).NE.ZERO) THEN
               IF(K.NE.NLAY .AND. LAYCBD(K).NE.0) THEN
                  IF(VKCB(J,I,LAYCBD(K)).NE.ZERO) GO TO 40
               END IF
               IF(K.NE.1 .AND. LAYCBD(K-1).NE.0) THEN
                  IF(VKCB(J,I,LAYCBD(K-1)).NE.ZERO) GO TO 40
               END IF
            END IF
         END IF
C
C2C-----ALL TRANSMISSIVE TERMS ARE ALL 0, SO CONVERT CELL TO NO FLOW.
         IBOUND(J,I,K)=0
         HNEW(J,I,K)=HCNV
         WETDRY(J,I,LAYWET(K))=ZERO
         WRITE(IOUT,43) K,I,J
   40    CONTINUE
C
      ELSE
C
C3------WETTING IS INACTIVE
         DO 50 I=1,NROW
         DO 50 J=1,NCOL
         IF(IBOUND(J,I,K).EQ.0) GO TO 50
C
C3A-----CHECK HORIZONTAL HYDRAULIC CONDUCTIVITY (HK).
         IF(HK(J,I,K).NE.ZERO) GO TO 50
C
C3B-----CHECK VERTICAL HYDRAULIC CONDUCTIVITY AND CONFINING BED
C3B-----VERTICAL HYDRAULIC CONDUCTIVITY.
         IF(NLAY.GT.1) THEN
            IF(VKA(J,I,K).NE.ZERO) THEN
               IF(K.NE.NLAY .AND. LAYCBD(K).NE.0) THEN
                  IF(VKCB(J,I,LAYCBD(K)).NE.ZERO) GO TO 50
               END IF
               IF(K.NE.1 .AND. LAYCBD(K-1).NE.0) THEN
                  IF(VKCB(J,I,LAYCBD(K-1)).NE.ZERO) GO TO 50
               END IF
            END IF
         END IF
C
C3C-----ALL TRANSMISSIVE TERMS ARE 0, SO CONVERT CELL TO NO FLOW.
         IBOUND(J,I,K)=0
         HNEW(J,I,K)=HCNV
         WRITE(IOUT,43) K,I,J
   43    FORMAT(1X,'NODE (LAYER,ROW,COL) ',I3,2(1X,I5),
     1 ' ELIMINATED BECAUSE ALL HYDRAULIC',/,
     2 ' CONDUCTIVITIES TO NODE ARE 0')
   50    CONTINUE
      END IF
   60 CONTINUE
C
C4------CALCULATE HOR. CONDUCTANCE(CR AND CC) FOR CONSTANT T LAYERS.
      DO 70 K=1,NLAY
      KK=K
      IF(LAYTYP(K).EQ.0) THEN
         CALL SGWF1LPF1HCOND(HNEW,IBOUND,CR,CC,HK,HANI,DELR,DELC,BOTM,
     1      NBOTM,KK,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT,WETDRY,
     2      WETFCT,IWETIT,IHDWET,HDRY)
      END IF
   70 CONTINUE
C
C5------CALCULATE VERTICAL CONDUCTANCE BETWEEN CONFINED LAYERS.
      IF(NLAY.GT.1) THEN
         DO 10 K=1,NLAY-1
         KK=K
         IF(LAYTYP(K).EQ.0 .AND. LAYTYP(K+1).EQ.0)
     1      CALL SGWF1LPF1VCOND(CV,HK,VKA,VKCB,IBOUND,BOTM,NBOTM,KK,
     2         NCOL,NROW,NLAY,HNEW,DELR,DELC,IOUT)
   10    CONTINUE
      END IF
C
C7------RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1B(HNEW,IBOUND,CR,CC,CV,BOTM,NCOL,NROW,NLAY,
     1      KSTP,KPER,ILPFCB,BUFF,IOUT,ICLPFL,DELT,PERTIM,TOTIM,
     2      IDIR,IBDRET,ICHFLG,IC1,IC2,IR1,IR2,IL1,IL2,NBOTM)
C
C-----VERSION 11JAN2000 SGWF1LPF1B
C     ******************************************************************
C     COMPUTE FLOW BETWEEN ADJACENT CELLS IN A SUBREGION OF THE GRID
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*16 TEXT(3)
      DOUBLE PRECISION HNEW,HD
C
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),
     1     CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),
     2     CV(NCOL,NROW,NLAY), BOTM(NCOL,NROW,0:NBOTM),
     3     BUFF(NCOL,NROW,NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA TEXT(1),TEXT(2),TEXT(3)
     1 /'FLOW RIGHT FACE ','FLOW FRONT FACE ','FLOW LOWER FACE '/
C     ------------------------------------------------------------------
C
C1------IF CELL-BY-CELL FLOWS WILL BE SAVED IN A FILE, SET FLAG IBD.
C1------RETURN IF FLOWS ARE NOT BEING SAVED OR RETURNED.
      ZERO=0.
      IBD=0
      IF(ILPFCB.GT.0) IBD=ICLPFL
      IF(IBD.EQ.0 .AND. IBDRET.EQ.0) RETURN
C
C2------SET THE SUBREGION EQUAL TO THE ENTIRE GRID IF VALUES ARE BEING
C2------SAVED IN A FILE.
      IF(IBD.NE.0) THEN
         K1=1
         K2=NLAY
         I1=1
         I2=NROW
         J1=1
         J2=NCOL
      END IF
C
C3------TEST FOR DIRECTION OF CALCULATION;  IF NOT ACROSS COLUMNS, GO TO
C3------STEP 4.  IF ONLY 1 COLUMN, RETURN.
      IF(IDIR.NE.1) GO TO 405
      IF(NCOL.EQ.1) RETURN
C
C3A-----CALCULATE FLOW ACROSS COLUMNS (THROUGH RIGHT FACE).  IF NOT
C3A-----SAVING IN A FILE, SET THE SUBREGION.  CLEAR THE BUFFER.
      IF(IBD.EQ.0) THEN
         K1=IL1
         K2=IL2
         I1=IR1
         I2=IR2
         J1=IC1-1
         IF(J1.LT.1) J1=1
         J2=IC2
      END IF
      DO 310 K=K1,K2
      DO 310 I=I1,I2
      DO 310 J=J1,J2
      BUFF(J,I,K)=ZERO
  310 CONTINUE
C
C3B-----FOR EACH CELL CALCULATE FLOW THRU RIGHT FACE & STORE IN BUFFER.
      IF(J2.EQ.NCOL) J2=J2-1
      DO 400 K=K1,K2
      DO 400 I=I1,I2
      DO 400 J=J1,J2
      IF(ICHFLG.EQ.0) THEN
         IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J+1,I,K).LE.0)) GO TO 400
      ELSE
         IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J+1,I,K).EQ.0)) GO TO 400
      END IF
      HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
      BUFF(J,I,K)=HDIFF*CR(J,I,K)
  400 CONTINUE
C
C3C-----RECORD CONTENTS OF BUFFER AND RETURN.
      IF(IBD.EQ.1)
     1   CALL UBUDSV(KSTP,KPER,TEXT(1),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(1),ILPFCB,BUFF,NCOL,NROW,
     1     NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      RETURN
C
C4------TEST FOR DIRECTION OF CALCULATION;  IF NOT ACROSS ROWS, GO TO
C4------STEP 5.  IF ONLY 1 ROW, RETURN.
  405 IF(IDIR.NE.2) GO TO 505
      IF(NROW.EQ.1) RETURN
C
C4A-----CALCULATE FLOW ACROSS ROWS (THROUGH FRONT FACE).  IF NOT SAVING
C4A-----IN A FILE, SET THE SUBREGION.  CLEAR THE BUFFER.
      IF(IBD.EQ.0) THEN
         K1=IL1
         K2=IL2
         I1=IR1-1
         IF(I1.LT.1) I1=1
         I2=IR2
         J1=IC1
         J2=IC2
      END IF
      DO 410 K=K1,K2
      DO 410 I=I1,I2
      DO 410 J=J1,J2
      BUFF(J,I,K)=ZERO
  410 CONTINUE
C
C4B-----FOR EACH CELL CALCULATE FLOW THRU FRONT FACE & STORE IN BUFFER.
      IF(I2.EQ.NROW) I2=I2-1
      DO 500 K=K1,K2
      DO 500 I=I1,I2
      DO 500 J=J1,J2
      IF(ICHFLG.EQ.0) THEN
         IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I+1,K).LE.0)) GO TO 500
      ELSE
         IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J,I+1,K).EQ.0)) GO TO 500
      END IF
      HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
      BUFF(J,I,K)=HDIFF*CC(J,I,K)
  500 CONTINUE
C
C4C-----RECORD CONTENTS OF BUFFER AND RETURN.
      IF(IBD.EQ.1)
     1   CALL UBUDSV(KSTP,KPER,TEXT(2),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(2),ILPFCB,BUFF,NCOL,NROW,
     1     NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      RETURN
C
C5------DIRECTION OF CALCULATION IS ACROSS LAYERS BY ELIMINATION.  IF
C5------ONLY 1 LAYER, RETURN.
  505 IF(NLAY.EQ.1) RETURN
C
C5A-----CALCULATE FLOW ACROSS LAYERS (THROUGH LOWER FACE).  IF NOT
C5A-----SAVING IN A FILE, SET THE SUBREGION.  CLEAR THE BUFFER.
      IF(IBD.EQ.0) THEN
         K1=IL1-1
         IF(K1.LT.1) K1=1
         K2=IL2
         I1=IR1
         I2=IR2
         J1=IC1
         J2=IC2
      END IF
      DO 510 K=K1,K2
      DO 510 I=I1,I2
      DO 510 J=J1,J2
      BUFF(J,I,K)=ZERO
  510 CONTINUE
C
C5B-----FOR EACH CELL CALCULATE FLOW THRU LOWER FACE & STORE IN BUFFER.
      IF(K2.EQ.NLAY) K2=K2-1
      DO 600 K=1,K2
      IF(K.LT.K1) GO TO 600
      DO 590 I=I1,I2
      DO 590 J=J1,J2
      IF(ICHFLG.EQ.0) THEN
         IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I,K+1).LE.0)) GO TO 590
      ELSE
         IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J,I,K+1).EQ.0)) GO TO 590
      END IF
      HD=HNEW(J,I,K+1)
      IF(LAYTYP(K+1).EQ.0) GO TO 580
      TMP=HD
      TOP=BOTM(J,I,LBOTM(K+1)-1)
      IF(TMP.LT.TOP) HD=TOP
  580 HDIFF=HNEW(J,I,K)-HD
      BUFF(J,I,K)=HDIFF*CV(J,I,K)
  590 CONTINUE
  600 CONTINUE
C
C5C-----RECORD CONTENTS OF BUFFER AND RETURN.
      IF(IBD.EQ.1)
     1   CALL UBUDSV(KSTP,KPER,TEXT(3),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(3),ILPFCB,BUFF,NCOL,NROW,
     1     NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      RETURN
      END
      SUBROUTINE SGWF1LPF1S(VBNM,VBVL,MSUM,HNEW,IBOUND,HOLD,SC1,
     1   BOTM,SC2,DELT,ISS,NCOL,NROW,NLAY,KSTP,KPER,ILPFCB,
     2   ICLPFL,BUFF,IOUT,PERTIM,TOTIM,NBOTM)
C-----VERSION 11JAN2000 SGWF1LPF1S
C     ******************************************************************
C     COMPUTE STORAGE BUDGET FLOW TERM FOR LPF.
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*16 VBNM(MSUM),TEXT
      DOUBLE PRECISION HNEW,STOIN,STOUT,SSTRG
C
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),
     1   HOLD(NCOL,NROW,NLAY),SC1(NCOL,NROW,NLAY),VBVL(4,MSUM),
     2   SC2(NCOL,NROW,NLAY),BOTM(NCOL,NROW,0:NBOTM),
     3   BUFF(NCOL,NROW,NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA TEXT /'         STORAGE'/
C     ------------------------------------------------------------------
C
C1------INITIALIZE BUDGET ACCUMULATORS AND 1/DELT.
      ZERO=0.
      STOIN=ZERO
      STOUT=ZERO
C2------IF STEADY STATE, STORAGE TERM IS ZERO
      IF(ISS.NE.0) GOTO 400
      ONE=1.
      TLED=ONE/DELT
C
C3------IF CELL-BY-CELL FLOWS WILL BE SAVED, SET FLAG IBD.
      IBD=0
      IF(ILPFCB.GT.0) IBD=ICLPFL
C
C4------CLEAR BUFFER.
      DO 210 K=1,NLAY
      DO 210 I=1,NROW
      DO 210 J=1,NCOL
      BUFF(J,I,K)=ZERO
210   CONTINUE
C
C5------LOOP THROUGH EVERY CELL IN THE GRID.
      KT=0
      DO 300 K=1,NLAY
      LC=LAYTYP(K)
      IF(LC.NE.0) KT=KT+1
      DO 300 I=1,NROW
      DO 300 J=1,NCOL
C
C6------SKIP NO-FLOW AND CONSTANT-HEAD CELLS.
      IF(IBOUND(J,I,K).LE.0) GO TO 300
      HSING=HNEW(J,I,K)
C
C7-----CHECK LAYER TYPE TO SEE IF ONE STORAGE CAPACITY OR TWO.
      IF(LC.EQ.0) GO TO 285
C
C7A----TWO STORAGE CAPACITIES.
      TP=BOTM(J,I,LBOTM(K)-1)
      RHO2=SC2(J,I,KT)*TLED
      RHO1=SC1(J,I,K)*TLED
      SOLD=RHO2
      IF(HOLD(J,I,K).GT.TP) SOLD=RHO1
      SNEW=RHO2
      IF(HSING.GT.TP) SNEW=RHO1
      STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING
      GO TO 288
C
C7B----ONE STORAGE CAPACITY.
  285 RHO=SC1(J,I,K)*TLED
      STRG=RHO*HOLD(J,I,K) - RHO*HSING
 
C
C8-----STORE CELL-BY-CELL FLOW IN BUFFER AND ADD TO ACCUMULATORS.
  288 BUFF(J,I,K)=STRG
      SSTRG=STRG
      IF(STRG) 292,300,294
  292 STOUT=STOUT-SSTRG
      GO TO 300
  294 STOIN=STOIN+SSTRG
C
  300 CONTINUE
C
C9-----IF IBD FLAG IS SET RECORD THE CONTENTS OF THE BUFFER.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,
     1                       ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT,ILPFCB,
     1            BUFF,NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
C
C10-----ADD TOTAL RATES AND VOLUMES TO VBVL & PUT TITLE IN VBNM.
  400 CONTINUE
      SIN=STOIN
      SOUT=STOUT
      VBVL(1,MSUM)=VBVL(1,MSUM)+SIN*DELT
      VBVL(2,MSUM)=VBVL(2,MSUM)+SOUT*DELT
      VBVL(3,MSUM)=SIN
      VBVL(4,MSUM)=SOUT
      VBNM(MSUM)=TEXT
      MSUM=MSUM+1
C
C11----RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1F(VBNM,VBVL,MSUM,HNEW,IBOUND,CR,CC,CV,BOTM,
     1         DELT,NCOL,NROW,NLAY,KSTP,KPER,ILPFCB,BUFF,IOUT,ICLPFL,
     2         PERTIM,TOTIM,NBOTM,ICHFLG)
C-----VERSION 13DEC2000 SGWF1LPF1F
C     ******************************************************************
C     COMPUTE FLOW FROM CONSTANT-HEAD CELLS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*16 VBNM(MSUM),TEXT
      DOUBLE PRECISION HNEW,HD,CHIN,CHOUT,XX1,XX2,XX3,XX4,XX5,XX6
C
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),
     1     CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),
     2     CV(NCOL,NROW,NLAY), VBVL(4,MSUM),
     3     BOTM(NCOL,NROW,0:NBOTM),BUFF(NCOL,NROW,NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
      DATA TEXT /'   CONSTANT HEAD'/
C     ------------------------------------------------------------------
C
C1------SET IBD TO INDICATE IF CELL-BY-CELL BUDGET VALUES WILL BE SAVED.
      IBD=0
      IF(ILPFCB.LT.0 .AND. ICLPFL.NE.0) IBD=-1
      IF(ILPFCB.GT.0) IBD=ICLPFL
C
C2------CLEAR BUDGET ACCUMULATORS.
      ZERO=0.
      CHIN=ZERO
      CHOUT=ZERO
      IBDLBL=0
C
C3------CLEAR BUFFER.
      DO 5 K=1,NLAY
      DO 5 I=1,NROW
      DO 5 J=1,NCOL
      BUFF(J,I,K)=ZERO
5     CONTINUE
C
C3A-----IF SAVING CELL-BY-CELL FLOW IN A LIST, COUNT CONSTANT-HEAD
C3A-----CELLS AND WRITE HEADER RECORDS.
      IF(IBD.EQ.2) THEN
         NCH=0
         DO 7 K=1,NLAY
         DO 7 I=1,NROW
         DO 7 J=1,NCOL
         IF(IBOUND(J,I,K).LT.0) NCH=NCH+1
7        CONTINUE
         CALL UBDSV2(KSTP,KPER,TEXT,ILPFCB,NCOL,NROW,NLAY,
     1          NCH,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      END IF
C
C4------LOOP THROUGH EACH CELL AND CALCULATE FLOW INTO MODEL FROM EACH
C4------CONSTANT-HEAD CELL.
      DO 200 K=1,NLAY
      DO 200 I=1,NROW
      DO 200 J=1,NCOL
C
C5------IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.
      IF (IBOUND(J,I,K).GE.0)GO TO 200
C
C6------CLEAR VALUES FOR FLOW RATE THROUGH EACH FACE OF CELL.
      X1=ZERO
      X2=ZERO
      X3=ZERO
      X4=ZERO
      X5=ZERO
      X6=ZERO
      CHCH1=ZERO
      CHCH2=ZERO
      CHCH3=ZERO
      CHCH4=ZERO
      CHCH5=ZERO
      CHCH6=ZERO
C
C7------CALCULATE FLOW THROUGH THE LEFT FACE.
C7------COMMENTS A-C APPEAR ONLY IN THE SECTION HEADED BY COMMENT 7,
C7------BUT THEY APPLY IN A SIMILAR MANNER TO SECTIONS 8-12.
C
C7A-----IF THERE IS NO FLOW TO CALCULATE THROUGH THIS FACE, THEN GO ON
C7A-----TO NEXT FACE.  NO FLOW OCCURS AT THE EDGE OF THE GRID, TO AN
C7A-----ADJACENT NO-FLOW CELL, OR TO AN ADJACENT CONSTANT-HEAD CELL.
      IF(J.EQ.1) GO TO 30
      IF(IBOUND(J-1,I,K).EQ.0) GO TO 30
      IF(IBOUND(J-1,I,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 30
C
C7B-----CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.
      HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)
      CHCH1=HDIFF*CR(J-1,I,K)
      IF(IBOUND(J-1,I,K).LT.0) GO TO 30
      X1=CHCH1
      XX1=X1
C
C7C-----ACCUMULATE POSITIVE AND NEGATIVE FLOW.
      IF (X1) 10,30,20
   10 CHOUT=CHOUT-XX1
      GO TO 30
   20 CHIN=CHIN+XX1
C
C8------CALCULATE FLOW THROUGH THE RIGHT FACE.
   30 IF(J.EQ.NCOL) GO TO 60
      IF(IBOUND(J+1,I,K).EQ.0) GO TO 60
      IF(IBOUND(J+1,I,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 60
      HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
      CHCH2=HDIFF*CR(J,I,K)
      IF(IBOUND(J+1,I,K).LT.0) GO TO 60
      X2=CHCH2
      XX2=X2
      IF(X2)40,60,50
   40 CHOUT=CHOUT-XX2
      GO TO 60
   50 CHIN=CHIN+XX2
C
C9------CALCULATE FLOW THROUGH THE BACK FACE.
   60 IF(I.EQ.1) GO TO 90
      IF (IBOUND(J,I-1,K).EQ.0) GO TO 90
      IF (IBOUND(J,I-1,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 90
      HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)
      CHCH3=HDIFF*CC(J,I-1,K)
      IF(IBOUND(J,I-1,K).LT.0) GO TO 90
      X3=CHCH3
      XX3=X3
      IF(X3) 70,90,80
   70 CHOUT=CHOUT-XX3
      GO TO 90
   80 CHIN=CHIN+XX3
C
C10-----CALCULATE FLOW THROUGH THE FRONT FACE.
   90 IF(I.EQ.NROW) GO TO 120
      IF(IBOUND(J,I+1,K).EQ.0) GO TO 120
      IF(IBOUND(J,I+1,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 120
      HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
      CHCH4=HDIFF*CC(J,I,K)
      IF(IBOUND(J,I+1,K).LT.0) GO TO 120
      X4=CHCH4
      XX4=X4
      IF (X4) 100,120,110
  100 CHOUT=CHOUT-XX4
      GO TO 120
  110 CHIN=CHIN+XX4
C
C11-----CALCULATE FLOW THROUGH THE UPPER FACE.
  120 IF(K.EQ.1) GO TO 150
      IF (IBOUND(J,I,K-1).EQ.0) GO TO 150
      IF (IBOUND(J,I,K-1).LT.0 .AND. ICHFLG.EQ.0) GO TO 150
      HD=HNEW(J,I,K)
      IF(LAYTYP(K).EQ.0) GO TO 122
      TMP=HD
      TOP=BOTM(J,I,LBOTM(K)-1)
      IF(TMP.LT.TOP) HD=TOP
  122 HDIFF=HD-HNEW(J,I,K-1)
      CHCH5=HDIFF*CV(J,I,K-1)
      IF(IBOUND(J,I,K-1).LT.0) GO TO 150
      X5=CHCH5
      XX5=X5
      IF(X5) 130,150,140
  130 CHOUT=CHOUT-XX5
      GO TO 150
  140 CHIN=CHIN+XX5
C
C12-----CALCULATE FLOW THROUGH THE LOWER FACE.
  150 IF(K.EQ.NLAY) GO TO 180
      IF(IBOUND(J,I,K+1).EQ.0) GO TO 180
      IF(IBOUND(J,I,K+1).LT.0 .AND. ICHFLG.EQ.0) GO TO 180
      HD=HNEW(J,I,K+1)
      IF(LAYTYP(K+1).EQ.0) GO TO 152
      TMP=HD
      TOP=BOTM(J,I,LBOTM(K+1)-1)
      IF(TMP.LT.TOP) HD=TOP
  152 HDIFF=HNEW(J,I,K)-HD
      CHCH6=HDIFF*CV(J,I,K)
      IF(IBOUND(J,I,K+1).LT.0) GO TO 180
      X6=CHCH6
      XX6=X6
      IF(X6) 160,180,170
  160 CHOUT=CHOUT-XX6
      GO TO 180
  170 CHIN=CHIN+XX6
C
C13-----SUM THE FLOWS THROUGH SIX FACES OF CONSTANT HEAD CELL, AND
C13-----STORE SUM IN BUFFER.
 180  RATE=CHCH1+CHCH2+CHCH3+CHCH4+CHCH5+CHCH6
      BUFF(J,I,K)=RATE
C
C14-----PRINT THE FLOW FOR THE CELL IF REQUESTED.
      IF(IBD.LT.0) THEN
         IF(IBDLBL.EQ.0) WRITE(IOUT,899) TEXT,KPER,KSTP
  899    FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
         WRITE(IOUT,900) K,I,J,RATE
  900    FORMAT(1X,'LAYER ',I3,'   ROW ',I5,'   COL ',I5,
     1       '   RATE ',1PG15.6)
         IBDLBL=1
      END IF
C
C15-----IF SAVING CELL-BY-CELL FLOW IN LIST, WRITE FLOW FOR CELL.
      IF(IBD.EQ.2) CALL UBDSVA(ILPFCB,NCOL,NROW,J,I,K,RATE,IBOUND,NLAY)
  200 CONTINUE
C
C16-----IF SAVING CELL-BY-CELL FLOW IN 3-D ARRAY, WRITE THE ARRAY.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,
     1                   ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT)
C
C17-----SAVE TOTAL CONSTANT HEAD FLOWS AND VOLUMES IN VBVL TABLE
C17-----FOR INCLUSION IN BUDGET. PUT LABELS IN VBNM TABLE.
      CIN=CHIN
      COUT=CHOUT
      VBVL(1,MSUM)=VBVL(1,MSUM)+CIN*DELT
      VBVL(2,MSUM)=VBVL(2,MSUM)+COUT*DELT
      VBVL(3,MSUM)=CIN
      VBVL(4,MSUM)=COUT
      VBNM(MSUM)=TEXT
      MSUM=MSUM+1
C
C18-----RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1SC(SC,NCOL,NROW,K,BOTM,NBOTM,DELR,DELC,ISPST)
C
C-----VERSION 11JAN2000 SGWF1LPF1SC
C     ******************************************************************
C     COMPUTE STORAGE CAPACITY
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION SC(NCOL,NROW),BOTM(NCOL,NROW,0:NBOTM),
     1          DELR(NCOL),DELC(NROW)
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
C     ------------------------------------------------------------------
C
C6A-----MULTIPLY SPECIFIC STORAGE BY THICKNESS, DELR, AND DELC TO GET
C6A-----CONFINED STORAGE CAPACITY.
      IF(ISPST.NE.0) THEN
         DO 80 I=1,NROW
         DO 80 J=1,NCOL
         THICK=BOTM(J,I,LBOTM(K)-1)-BOTM(J,I,LBOTM(K))
         SC(J,I)=SC(J,I)*THICK*DELR(J)*DELC(I)
   80    CONTINUE
      ELSE
C
C6B-----IF LAYER IS CONVERTIBLE MULTIPLY SPECIFIC YIELD BY DELR AND DELC
C6B-----TO GET UNCONFINED STORAGE CAPACITY(SC2).
         DO 85 I=1,NROW
         DO 85 J=1,NCOL
         SC(J,I)=SC(J,I)*DELR(J)*DELC(I)
   85    CONTINUE
      END IF
C
      RETURN
      END
      SUBROUTINE SGWF1LPF1HCOND(HNEW,IBOUND,CR,CC,HK,HANI,DELR,DELC,
     1 BOTM,NBOTM,K,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT,WETDRY,
     2 WETFCT,IWETIT,IHDWET,HDRY)
C-----VERSION 13JAN2000 SGWF1LPF1HCOND
C     ******************************************************************
C     COMPUTE HORIZONTAL BRANCH CONDUCTANCE FOR ONE LAYER.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------

      use altparam                               !jd

      DOUBLE PRECISION HNEW,HHD,BBOT,TTOP
C
      DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     1 CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY), HK(NCOL,NROW,NLAY),
     2 HANI(NCOL,NROW,NLAY), DELR(NCOL), DELC(NROW),
     3 BOTM(NCOL,NROW,0:NBOTM),WETDRY(NCOL,NROW,NLAY)
      CHARACTER*3 ACNVRT
      DIMENSION ICNVRT(5),JCNVRT(5),ACNVRT(5)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
C1------INITIALIZE DATA.
      ZERO=0.
      NCNVRT=0
      IHDCNV=0
C
C2------IF LAYER IS WETTABLE CONVERT DRY CELLS TO WET WHEN APPROPRIATE.
      ITFLG=1
      IF(LAYWET(K).NE.0) ITFLG=MOD(KITER,IWETIT)
      IF(ITFLG.EQ.0) CALL SGWF1LPF1WET(HNEW,IBOUND,BOTM,NBOTM,K,KITER,
     1      KSTP,KPER,NCOL,NROW,NLAY,IOUT,WETDRY,WETFCT,IHDWET,
     2      IHDCNV,NCNVRT,ICNVRT,JCNVRT,ACNVRT)
C
C3------LOOP THROUGH EACH CELL, AND CALCULATE SATURATED THICKNESS.
      DO 200 I=1,NROW
      DO 200 J=1,NCOL
C
C3A-----SET STAURATED THICKNESS=0. FOR DRY CELLS.
      IF(IBOUND(J,I,K).EQ.0) THEN
         CC(J,I,K)=ZERO
      ELSE
C
C3B-----CALCULATE SATURATED THICKNESS FOR A WET CELL.
         BBOT=BOTM(J,I,LBOTM(K))
         TTOP=BOTM(J,I,LBOTM(K)-1)
         IF(BBOT.GT.TTOP) THEN
            WRITE(IOUT,35) K,I,J
   35       FORMAT(1X,/1X,'Negative cell thickness at (layer,row,col)',
     1      I4,',',I5,',',I5)
            WRITE(IOUT,36) TTOP,BBOT
   36       FORMAT(1X,'Top elevation, bottom elevation:',1P,2G13.5)
        call stopfile ! emrl
            STOP
         END IF
         IF(LAYTYP(K).NE.0) THEN
            HHD=HNEW(J,I,K)
            IF(HHD.LT.TTOP) TTOP=HHD
         END IF
         THCK=TTOP-BBOT
          if(minthick.gt.0.0)then                               !jd
           if((k.eq.nlay).and.(thck.le.minthick))thck=minthick !jd
         end if                                                !jd
         CC(J,I,K)=THCK
C
C
C3C-----WHEN SATURATED THICKNESS <= 0, PRINT A MESSAGE AND SET
C3C-----HNEW=HDRY, SATURATED THICKNESS=0.0, AND IBOUND=0.
         IF(THCK.LE.ZERO) THEN
            CALL SGWF1LPF1WDMSG(1,NCNVRT,ICNVRT,JCNVRT,ACNVRT,IHDCNV,
     1             IOUT,KITER,J,I,K,KSTP,KPER,NCOL,NROW)
            if(hdrybot.eq.0) then                              !jd
            HNEW(J,I,K)=HDRY
            else                                                !jd
            hnew(j,i,k)=bbot                                    !jd
            end if                                              !jd
            CC(J,I,K)=ZERO
            IF(IBOUND(J,I,K).LT.0) THEN
               WRITE(IOUT,151)
  151          FORMAT(1X,/1X,'CONSTANT-HEAD CELL WENT DRY',
     1          ' -- SIMULATION ABORTED')
                    WRITE(IOUT,*) TTOP, BBOT, THCK
               WRITE(IOUT,152) K,I,J,KITER,KSTP,KPER
  152          FORMAT(1X,'LAYER=',I3,' ROW=',I5,' COLUMN=',I5,
     1    ' ITERATION=',I3,' TIME STEP=',I3,' STRESS PERIOD=',I4)
        call stopfile ! emrl
               STOP
            END IF
            IBOUND(J,I,K)=0
         END IF
      END IF
  200 CONTINUE
C
C4------PRINT ANY REMAINING CELL CONVERSIONS NOT YET PRINTED.
      CALL SGWF1LPF1WDMSG(0,NCNVRT,ICNVRT,JCNVRT,ACNVRT,IHDCNV,
     1             IOUT,KITER,J,I,K,KSTP,KPER,NCOL,NROW)
C
C5------CHANGE IBOUND VALUE FOR CELLS THAT CONVERTED TO WET THIS
C5------ITERATION FROM 30000 to 1.
      IF(LAYWET(K).NE.0) THEN
         DO 205 I=1,NROW
         DO 205 J=1,NCOL
         IF(IBOUND(J,I,K).EQ.30000) IBOUND(J,I,K)=1
  205    CONTINUE
      END IF
C
C6------COMPUTE HORIZONTAL BRANCH CONDUCTANCES FROM CELL HYDRAULIC
C6------CONDUCTIVITY, SATURATED THICKNESS, AND GRID DIMENSIONS.
      IF(LAYAVG(K).EQ.0) THEN
         CALL SGWF1LPF1HHARM(CR,CC,HK,HANI,IBOUND,DELR,DELC,K,NCOL,NROW,
     1         NLAY)
      ELSE IF(LAYAVG(K).EQ.1) THEN
         CALL SGWF1LPF1HLOG(IBOUND,CR,CC,HK,HANI,DELR,DELC,K,NCOL,
     1        NROW,NLAY)
      ELSE
         CALL SGWF1LPF1HUNCNF(IBOUND,CR,CC,HK,HANI,DELR,DELC,K,NCOL,
     1        NROW,NLAY)
      END IF
C
C7------RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1WET(HNEW,IBOUND,BOTM,NBOTM,K,KITER,KSTP,KPER,
     1      NCOL,NROW,NLAY,IOUT,WETDRY,WETFCT,IHDWET,IHDCNV,
     2      NCNVRT,ICNVRT,JCNVRT,ACNVRT)
C
C-----VERSION 11JAN2000 SGWF1LPF1WET
C     ******************************************************************
C     CONVERT DRY CELLS TO WET.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW
C
      DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     1          BOTM(NCOL,NROW,0:NBOTM),WETDRY(NCOL,NROW,NLAY)
      CHARACTER*3 ACNVRT
      DIMENSION ICNVRT(5),JCNVRT(5),ACNVRT(5)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
      ZERO=0.0
C
C1------LOOP THROUGH ALL CELLS.
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
C
C2------IF CELL IS DRY AND IF IT IS WETTABLE, CONTINUE CHECKING TO SEE
C2------IF IT SHOULD BECOME WET.
      IF(IBOUND(J,I,K).EQ.0 .AND. WETDRY(J,I,LAYWET(K)).NE.ZERO) THEN
C
C3------CALCULATE WETTING ELEVATION.
         WD=WETDRY(J,I,LAYWET(K))
         IF(WD.LT.ZERO) WD=-WD
         TURNON=BOTM(J,I,LBOTM(K))+WD
C
C4------CHECK HEAD IN CELL BELOW TO SEE IF WETTING ELEVATION HAS BEEN
C4------REACHED.
         IF(K.NE.NLAY) THEN
            HTMP=HNEW(J,I,K+1)
            IF(IBOUND(J,I,K+1).GT.0 .AND. HTMP.GE.TURNON) GO TO 50
         END IF
C
C5------CHECK HEAD IN ADJACENT HORIZONTAL CELLS TO SEE IF WETTING
C5------ELEVATION HAS BEEN REACHED.
         IF(WETDRY(J,I,LAYWET(K)).GT.ZERO) THEN
            IF(J.NE.1) THEN
               HTMP=HNEW(J-1,I,K)
               IF(IBOUND(J-1,I,K).GT.0 .AND. IBOUND(J-1,I,K).NE.30000.
     1                       AND. HTMP.GE.TURNON) GO TO 50
            END IF
            IF(J.NE.NCOL) THEN
               HTMP=HNEW(J+1,I,K)
               IF(IBOUND(J+1,I,K).GT.0 .AND. HTMP.GE.TURNON) GO TO 50
            END IF
            IF(I.NE.1) THEN
               HTMP=HNEW(J,I-1,K)
               IF(IBOUND(J,I-1,K).GT.0 .AND. IBOUND(J,I-1,K).NE.30000.
     1                       AND. HTMP.GE.TURNON) GO TO 50
            END IF
            IF(I.NE.NROW) THEN
               HTMP=HNEW(J,I+1,K)
               IF(IBOUND(J,I+1,K).GT.0 .AND. HTMP.GE.TURNON) GO TO 50
            END IF
         END IF
C
C6------WETTING ELEVATION HAS NOT BEEN REACHED, SO CELL REMAINS DRY.
         GO TO 100
C
C7------CELL BECOMES WET.  PRINT MESSAGE, SET INITIAL HEAD, AND SET
C7------IBOUND.
   50    CALL SGWF1LPF1WDMSG(2,NCNVRT,ICNVRT,JCNVRT,ACNVRT,IHDCNV,
     1             IOUT,KITER,J,I,K,KSTP,KPER,NCOL,NROW)
C
C7A-----USE EQUATION 3A IF IHDWET=0; ISE EQUATION 3B IF IHDWET IS NOT 0.
         IF(IHDWET.EQ.0) THEN
            HNEW(J,I,K)=BOTM(J,I,LBOTM(K))+
     1                        WETFCT*(HTMP-BOTM(J,I,LBOTM(K)))
         ELSE
            HNEW(J,I,K)=BOTM(J,I,LBOTM(K))+WETFCT*WD
         END IF
         IBOUND(J,I,K)=30000
      END IF
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE SGWF1LPF1WDMSG(ICODE,NCNVRT,ICNVRT,JCNVRT,ACNVRT,
     1             IHDCNV,IOUT,KITER,J,I,K,KSTP,KPER,NCOL,NROW)
C
C-----VERSION 11JAN2000 SGWF1LPF1WDMSG
C     ******************************************************************
C     PRINT MESSAGE WHEN CELLS CONVERT BETWEEN WET AND DRY.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*3 ACNVRT
      DIMENSION ICNVRT(5),JCNVRT(5),ACNVRT(5)
C     ------------------------------------------------------------------
C1------KEEP TRACK OF CELL CONVERSIONS.
      IF(ICODE.GT.0) THEN
         NCNVRT=NCNVRT+1
         ICNVRT(NCNVRT)=I
         JCNVRT(NCNVRT)=J
         IF(ICODE.EQ.1) THEN
            ACNVRT(NCNVRT)='DRY'
         ELSE
            ACNVRT(NCNVRT)='WET'
         END IF
      END IF
C
C2------PRINT A LINE OF DATA IF 5 CONVERSIONS HAVE OCCURRED OR IF ICODE
C2------INDICATES THAT A PARTIAL LINE SHOULD BE PRINTED.
      IF(NCNVRT.EQ.5 .OR. (ICODE.EQ.0 .AND. NCNVRT.GT.0)) THEN
         IF(IHDCNV.EQ.0) WRITE(IOUT,17) KITER,K,KSTP,KPER
   17    FORMAT(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I3,'  LAYER=',
     1       I3,'  STEP=',I3,'  PERIOD=',I4,'   (ROW,COL)')
         IHDCNV=1
         IF (NROW.LE.999 .AND. NCOL.LE.999) THEN
           WRITE(IOUT,18) (ACNVRT(L),ICNVRT(L),JCNVRT(L),L=1,NCNVRT)
   18      FORMAT(1X,3X,5(A,'(',I3,',',I3,')   '))
         ELSE
           WRITE(IOUT,19) (ACNVRT(L),ICNVRT(L),JCNVRT(L),L=1,NCNVRT)
   19      FORMAT(1X,4(3X,A,'(',I5,',',I5,')'))
         ENDIF
         NCNVRT=0
      END IF
C
C3------RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1HHARM(CR,CC,HK,HANI,IBOUND,DELR,DELC,K,NCOL,
     1         NROW,NLAY)
C
C-----VERSION 11JAN2000 SGWF1LPF1HHARM
C     ******************************************************************
C     COMPUTE HORIZONTAL BRANCH CONDUCTANCE USING HARMONIC MEAN OF BLOCK
C     CONDUCTANCES (DISTANCE WEIGHTED HARMONIC MEAN OF TRANSMISSIVITY).
C     CELL THICKNESS IS IN CC UPON ENTRY.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     1          HK(NCOL,NROW,NLAY),HANI(NCOL,NROW,NLAY),
     2          IBOUND(NCOL,NROW,NLAY),DELR(NCOL),DELC(NROW)
C
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
C     ------------------------------------------------------------------
      ZERO=0.
      TWO=2.
C
C1------FOR EACH CELL CALCULATE BRANCH CONDUCTANCES FROM THAT CELL
C1------TO THE ONE ON THE RIGHT AND THE ONE IN FRONT.
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
C
C2------IF CELL IS DRY OR HK=0., SET CONDUCTANCE EQUAL TO 0 AND GO ON
C2------TO NEXT CELL.
      IF(IBOUND(J,I,K).EQ.0 .OR. HK(J,I,K).EQ.ZERO) THEN
         CR(J,I,K)=ZERO
         CC(J,I,K)=ZERO
      ELSE
C
C3------CELL IS WET -- CALCULATE TRANSMISSIVITY OF CELL.
         T1=HK(J,I,K)*CC(J,I,K)
C3A-----IF THIS IS NOT THE LAST COLUMN (RIGHTMOST), CALCULATE
C3A-----BRANCH CONDUCTANCE IN THE ROW DIRECTION (CR) TO THE RIGHT.
         IF(J.NE.NCOL) THEN
            IF(IBOUND(J+1,I,K).NE.0) THEN
               T2=HK(J+1,I,K)*CC(J+1,I,K)
               CR(J,I,K)=TWO*T2*T1*DELC(I)/(T1*DELR(J+1)+T2*DELR(J))
            ELSE
               CR(J,I,K)=ZERO
            END IF
         ELSE
C3B-----IF THIS IS THE LAST COLUMN, SET BRANCH CONDUCTANCE=0.
            CR(J,I,K)=ZERO
         END IF
C
C3C-----IF THIS IS NOT THE LAST ROW (FRONTMOST) THEN CALCULATE
C3C-----BRANCH CONDUCTANCE IN THE COLUMN DIRECTION (CC) TO THE FRONT.
         IF(I.NE.NROW) THEN
            IF(IBOUND(J,I+1,K).NE.0) THEN
               T2=HK(J,I+1,K)*CC(J,I+1,K)
               IF(CHANI(K).LE.ZERO) THEN
                  KHANI=-CHANI(K)
                  T1=T1*HANI(J,I,KHANI)
                  T2=T2*HANI(J,I+1,KHANI)
               ELSE
                  T1=T1*CHANI(K)
                  T2=T2*CHANI(K)
               END IF
               CC(J,I,K)=TWO*T2*T1*DELR(J)/(T1*DELC(I+1)+T2*DELC(I))
            ELSE
C3D-----IF THIS IS THE LAST ROW, SET BRANCH CONDUCTANCE=0.
               CC(J,I,K)=ZERO
            END IF
         ELSE
            CC(J,I,K)=ZERO
         END IF
      END IF
  100 CONTINUE
C
C4------RETURN
      RETURN
      END
      SUBROUTINE SGWF1LPF1HLOG(IBOUND,CR,CC,HK,HANI,DELR,DELC,K,NCOL,
     1       NROW,NLAY)
C
C-----VERSION 11JAN2000 SGWF1LPF1HLOG
C     ******************************************************************
C-----COMPUTE HORIZONTAL CONDUCTANCE USING LOGARITHMIC MEAN
C-----TRANSMISSIVITY -- ACTIVATED BY LAYAVG=1
C-----CELL SATURATED THICKNESS IS IN CC.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DIMENSION CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     1    HK(NCOL,NROW,NLAY),HANI(NCOL,NROW,NLAY),
     2    IBOUND(NCOL,NROW,NLAY),DELR(NCOL),DELC(NROW)
C
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
      ZERO=0.
      TWO=2.
      HALF=0.5
      FRAC1=1.005
      FRAC2=0.995
C
C1------FOR EACH CELL CALCULATE BRANCH CONDUCTANCES FROM THAT CELL
C1------TO THE ONE ON THE RIGHT AND THE ONE IN FRONT.
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
C
C2------IF CELL IS DRY OR HK=0., SET CONDUCTANCE EQUAL TO 0 AND GO ON
C2------TO NEXT CELL.
      IF(IBOUND(J,I,K).EQ.0 .OR. HK(J,I,K).EQ.ZERO) THEN
         CR(J,I,K)=ZERO
         CC(J,I,K)=ZERO
      ELSE
C
C3------CELL IS WET -- CALCULATE TRANSMISSIVITY OF CELL.
         T1=HK(J,I,K)*CC(J,I,K)
C3A-----IF THIS IS NOT THE LAST COLUMN(RIGHTMOST) THEN CALCULATE
C3A-----BRANCH CONDUCTANCE IN THE ROW DIRECTION (CR) TO THE RIGHT.
         IF(J.NE.NCOL) THEN
            IF(IBOUND(J+1,I,K).NE.0) THEN
C3A1----LOGARITHMIC MEAN INTERBLOCK TRANSMISSIVITY
               T2=HK(J+1,I,K)*CC(J+1,I,K)
               RATIO=T2/T1
               IF(RATIO.GT.FRAC1 .OR. RATIO.LT.FRAC2) THEN
                  T=(T2-T1)/LOG(RATIO)
               ELSE
                  T=HALF*(T1+T2)
               END IF
               CR(J,I,K)=TWO*DELC(I)*T/(DELR(J+1)+DELR(J))
            ELSE
               CR(J,I,K)=ZERO
            END IF
         ELSE
            CR(J,I,K)=ZERO
         END IF
C
C3B-----IF THIS IS NOT THE LAST ROW (FRONTMOST) THEN CALCULATE
C3B-----BRANCH CONDUCTANCE IN THE COLUMN DIRECTION (CC) TO THE FRONT.
         IF(I.NE.NROW) THEN
            IF(IBOUND(J,I+1,K).NE.0) THEN
               T2=HK(J,I+1,K)*CC(J,I+1,K)
               IF(CHANI(K).LE.ZERO) THEN
                  KHANI=-CHANI(K)
                  T1=T1*HANI(J,I,KHANI)
                  T2=T2*HANI(J,I+1,KHANI)
               ELSE
                  T1=T1*CHANI(K)
                  T2=T2*CHANI(K)
               END IF
               RATIO=T2/T1
               IF(RATIO.GT.FRAC1 .OR. RATIO.LT.FRAC2) THEN
                  T=(T2-T1)/LOG(RATIO)
               ELSE
                  T=HALF*(T1+T2)
               END IF
               CC(J,I,K)=TWO*DELR(J)*T/(DELC(I+1)+DELC(I))
            ELSE
               CC(J,I,K)=ZERO
            END IF
         ELSE
            CC(J,I,K)=ZERO
         END IF
      END IF
  100 CONTINUE
C
C4------RETURN
      RETURN
      END
      SUBROUTINE SGWF1LPF1HUNCNF(IBOUND,CR,CC,HK,HANI,DELR,DELC,K,NCOL,
     1        NROW,NLAY)
C
C-----VERSION 11JAN2000 SGWF1LPF1HUNCNF
C     ******************************************************************
C-----COMPUTE HORIZONTAL CONDUCTANCE USING ARITHMETIC MEAN SATURATED
C-----THICKNESS AND LOGARITHMIC MEAN HYDRAULIC CONDUCTIVITY.
C-----CELL SATURATED THICKNESS IS IN CC.
C-----ACTIVATED BY LAYAVG=2
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DIMENSION CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     1   HK(NCOL,NROW,NLAY),HANI(NCOL,NROW,NLAY),
     2   IBOUND(NCOL,NROW,NLAY),DELR(NCOL),DELC(NROW)
C
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
      ZERO=0.
      HALF=0.5
      FRAC1=1.005
      FRAC2=0.995
C
C1------FOR EACH CELL CALCULATE BRANCH CONDUCTANCES FROM THAT CELL
C1------TO THE ONE ON THE RIGHT AND THE ONE IN FRONT.
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
C
C2------IF CELL IS DRY OR HK=0., SET CONDUCTANCE EQUAL TO 0 AND GO ON
C2------TO NEXT CELL.
      IF(IBOUND(J,I,K).EQ.0 .OR. HK(J,I,K).EQ.ZERO) THEN
         CR(J,I,K)=ZERO
         CC(J,I,K)=ZERO
      ELSE
C
C3------CELL IS WET -- CALCULATE TRANSMISSIVITY OF CELL.
         HYC1=HK(J,I,K)
C3A-----IF THIS IS NOT THE LAST COLUMN(RIGHTMOST) THEN CALCULATE
C3A-----BRANCH CONDUCTANCE IN THE ROW DIRECTION (CR) TO THE RIGHT.
         IF(J.NE.NCOL) THEN
            IF(IBOUND(J+1,I,K).NE.0) THEN
C3A1----LOGARITHMIC MEAN HYDRAULIC CONDUCTIVITY
               HYC2=HK(J+1,I,K)
               RATIO=HYC2/HYC1
               IF(RATIO.GT.FRAC1 .OR. RATIO.LT.FRAC2) THEN
                  HYC=(HYC2-HYC1)/LOG(RATIO)
               ELSE
                  HYC=HALF*(HYC1+HYC2)
               END IF
C3A2----MULTIPLY LOGARITHMIC K BY ARITMETIC SATURATED THICKNESS.
               CR(J,I,K)=DELC(I)*HYC*(CC(J,I,K)+CC(J+1,I,K))/
     1                      (DELR(J+1)+DELR(J))
            ELSE
               CR(J,I,K)=ZERO
            END IF
         ELSE
            CR(J,I,K)=ZERO
         END IF
C
C3B-----IF THIS IS NOT THE LAST ROW (FRONTMOST) THEN CALCULATE
C3B-----BRANCH CONDUCTANCE IN THE COLUMN DIRECTION (CC) TO THE FRONT.
         IF(I.NE.NROW) THEN
            IF(IBOUND(J,I+1,K).NE.0) THEN
C3B1----LOGARITHMIC MEAN HYDRAULIC CONDUCTIVITY
               HYC2=HK(J,I+1,K)
               IF(CHANI(K).LE.ZERO) THEN
                  KHANI=-CHANI(K)
                  HYC1=HYC1*HANI(J,I,KHANI)
                  HYC2=HYC2*HANI(J,I+1,KHANI)
               ELSE
                  HYC1=HYC1*CHANI(K)
                  HYC2=HYC2*CHANI(K)
               END IF
               RATIO=HYC2/HYC1
               IF(RATIO.GT.FRAC1 .OR. RATIO.LT.FRAC2) THEN
                  HYC=(HYC2-HYC1)/LOG(RATIO)
               ELSE
                  HYC=HALF*(HYC1+HYC2)
               END IF
C3B2----MULTIPLY LOGARITHMIC K BY ARITMETIC SATURATED THICKNESS.
               CC(J,I,K)=DELR(J)*HYC*(CC(J,I,K)+CC(J,I+1,K))/
     1                      (DELC(I+1)+DELC(I))
            ELSE
               CC(J,I,K)=ZERO
            END IF
         ELSE
            CC(J,I,K)=ZERO
         END IF
      END IF
  100 CONTINUE
C
C4------RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1VCOND(CV,HK,VKA,VKCB,IBOUND,BOTM,NBOTM,K,
     1         NCOL,NROW,NLAY,HNEW,DELR,DELC,IOUT)
C
C-----VERSION 13JAN2000 SGWF1LPF1VHARM
C     ******************************************************************
C     COMPUTE VERTICAL BRANCH CONDUCTANCE BETWEEN A LAYER AND THE NEXT
C     LOWER LAYER FROM VERTICAL HYDRAULIC CONDUCTIVITY.
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW,BBOT,TTOP,HHD
C
      DIMENSION CV(NCOL,NROW,NLAY),HK(NCOL,NROW,NLAY),
     1      VKA(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     2      BOTM(NCOL,NROW,0:NBOTM),HNEW(NCOL,NROW,NLAY),
     3      VKCB(NCOL,NROW,NLAY),DELR(NCOL),DELC(NROW)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C
C     ------------------------------------------------------------------
      IF(K.EQ.NLAY) RETURN
      ZERO=0.
      HALF=0.5
C
C1------LOOP THROUGH ALL CELLS IN THE LAYER.
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
      CV(J,I,K)=ZERO
      IF(IBOUND(J,I,K).NE.0 .AND. IBOUND(J,I,K+1).NE.0) THEN
C
C2------CALCULATE VERTICAL HYDRAULIC CONDUCTIVITY FOR CELL.
         IF(LAYVKA(K).EQ.0) THEN
            HYC1=VKA(J,I,K)
         ELSE
            HYC1=HK(J,I,K)/VKA(J,I,K)
         END IF
         IF(HYC1.GT.ZERO) THEN
C3------CALCULATE VERTICAL HYDRAULIC CONDUCTIVITY FOR CELL BELOW.
            IF(LAYVKA(K+1).EQ.0) THEN
               HYC2=VKA(J,I,K+1)
            ELSE
               HYC2=(HK(J,I,K+1)/VKA(J,I,K+1))
            END IF
            IF(HYC2.GT.ZERO) THEN
C
C4------CALCULATE INVERSE LEAKANCE FOR CELL.
               BBOT=BOTM(J,I,LBOTM(K))
               TTOP=BOTM(J,I,LBOTM(K)-1)
               IF(LAYTYP(K).NE.0) THEN
                  HHD=HNEW(J,I,K)
                  IF(HHD.LT.TTOP) TTOP=HHD
               END IF
               BOVK1=(TTOP-BBOT)*HALF/HYC1
C
C5------CALCULATE INVERSE LEAKANCE FOR CELL BELOW.
               BBOT=BOTM(J,I,LBOTM(K+1))
               TTOP=BOTM(J,I,LBOTM(K+1)-1)
               B=(TTOP-BBOT)*HALF
               IF(LAYTYP(K+1).NE.0) THEN
                  HHD=HNEW(J,I,K+1)
                  IF(HHD.LT.TTOP) B=ZERO
               END IF
               BOVK2=B/HYC2
C
               IF(LAYCBD(K).NE.0) THEN
C6------CALCULATE VERTICAL HYDRAULIC CONDUCTIVITY FOR CONFINING BED.
                  IF(VKCB(J,I,LAYCBD(K)).GT.ZERO) THEN
C
C7------CALCULATE INVERSE LEAKANCE FOR CONFINING BED.
                     B=BOTM(J,I,LBOTM(K))-BOTM(J,I,LBOTM(K)+1)
                     IF(B.LT.ZERO) THEN
                        WRITE(IOUT,45) K,I,J
   45                   FORMAT(1X,/1X,
     1  'Negative confining bed thickness below cell (Layer,row,col)',
     2                  I4,',',I5,',',I5)
            WRITE(IOUT,46) BOTM(J,I,LBOTM(K)),BOTM(J,I,LBOTM(K)+1)
   46       FORMAT(1X,'Top elevation, bottom elevation:',1P,2G13.5)
        call stopfile ! emrl
                        STOP
                     END IF
                     CBBOVK=B/VKCB(J,I,LAYCBD(K))
                     CV(J,I,K)=DELR(J)*DELC(I)/(BOVK1+CBBOVK+BOVK2)
                  END IF
               ELSE
                  CV(J,I,K)=DELR(J)*DELC(I)/(BOVK1+BOVK2)
               END IF
            END IF
         END IF
      END IF
  100 CONTINUE
C
C8------RETURN.
      RETURN
      END
      SUBROUTINE SGWF1LPF1CK(IOUT,NP,PTYP)
C-----VERSION 11JAN2000 SGWF1LPF1CK
C     ******************************************************************
C     CHECK THAT JUST-DEFINED PARAMETER OF TYPE 'VK' OR 'VANI' IS USED
C     CONSISTENTLY WITH LAYVKA ENTRIES FOR LAYERS LISTED IN CLUSTERS FOR
C     THE PARAMETER
C     ******************************************************************
C
C      SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*4 PTYP
      INCLUDE 'param.inc'
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
  590 FORMAT(/,
     &1X,'LAYVKA entered for layer ',i3,' is: ',i3,'; however,',
     &' layer ',i3,' is',/,' listed in a cluster for parameter "',a,
     &'" of type ',a,' and')
  600 FORMAT(
     &1X,'parameters of type VK can apply only to layers for which',
     &/,' LAYVKA is specified as zero -- STOP EXECUTION (SGWF1LPF1CK)')
  610 FORMAT(
     &1X,'parameters of type VANI can apply only to layers for which',/,
     &' LAYVKA is not specified as zero -- STOP EXECUTION',
     &' (SGWF1LPF1CK)')
C
C     LOOP THROUGH CLUSTERS FOR THIS PARAMETER
      DO 10 ICL = IPLOC(1,NP),IPLOC(2,NP)
        LAY = IPCLST(1,ICL)
        LV = LAYVKA(LAY)
        IF (PTYP.EQ.'VK  ' .AND. LV.NE.0) THEN
          WRITE (IOUT,590) LAY,LV,LAY,PARNAM(NP),'VK'
          WRITE (IOUT,600)
        call stopfile ! emrl
          STOP
        ELSEIF (PTYP.EQ.'VANI' .AND. LV.EQ.0) THEN
          WRITE (IOUT,590) LAY,LV,LAY,PARNAM(NP),'VANI'
          WRITE (IOUT,610)
        call stopfile ! emrl
          STOP
        ENDIF
   10 CONTINUE
C
      RETURN
      END
